Precipation-2D

import intake
from easygems import healpix as egh
import easygems.remap as egr
import healpy as hp

import matplotlib.pyplot as plt
import warnings

warnings.filterwarnings("ignore", category=FutureWarning) # don't warn us about future package conflicts

open dataset

current_location = "NCAR" #"online"
cat = intake.open_catalog("https://digital-earths-global-hackathon.github.io/catalog/catalog.yaml")[current_location]
list(cat)
['CERES_EBAF',
 'ERA5',
 'IR_IMERG',
 'JRA3Q',
 'MERRA2',
 'casesm2_10km_nocumulus',
 'icon_d3hp003',
 'icon_d3hp003aug',
 'icon_d3hp003feb',
 'icon_ngc4008',
 'ifs_tco3999-ng5_deepoff',
 'ifs_tco3999-ng5_rcbmf',
 'ifs_tco3999-ng5_rcbmf_cf',
 'mpas_dyamond1',
 'mpas_dyamond2',
 'mpas_dyamond3',
 'nicam_220m_test',
 'nicam_gl11',
 'scream-dkrz',
 'scream2D_hrly',
 'scream_lnd',
 'scream_ne120',
 'tracking-d3hp003',
 'um_Africa_km4p4_RAL3P3_n1280_GAL9_nest',
 'um_CTC_km4p4_RAL3P3_n1280_GAL9_nest',
 'um_SAmer_km4p4_RAL3P3_n1280_GAL9_nest',
 'um_SEA_km4p4_RAL3P3_n1280_GAL9_nest',
 'um_glm_n1280_CoMA9_TBv1p2',
 'um_glm_n1280_GAL9',
 'um_glm_n2560_RAL3p3',
 'wrf_conus',
 'wrf_samerica']
import pandas as pd

pd.DataFrame(cat["mpas_dyamond1"].describe()["user_parameters"])
name description type allowed default
0 time temporal resolution of the dataset str [PT15M, PT3H] PT3H
1 zoom zoom resolution of the dataset int [10, 9, 8, 7, 6, 5, 4, 3, 2, 1] 7

select variable

scream = cat['scream2D_hrly'](zoom = 7).to_dask()
scream
<xarray.Dataset> Size: 15GB
Dimensions:   (crs: 0, time: 9457, cell: 196608)
Coordinates:
  * crs       (crs) float64 0B 
  * time      (time) object 76kB 2019-08-01 00:00:00 ... 2020-09-01 00:00:00
Dimensions without coordinates: cell
Data variables:
    pr        (time, cell) float32 7GB dask.array<chunksize=(24, 196608), meta=np.ndarray>
    rlut      (time, cell) float32 7GB dask.array<chunksize=(24, 196608), meta=np.ndarray>
    ELEV      (cell) float32 786kB dask.array<chunksize=(4096,), meta=np.ndarray>
    LANDFRAC  (cell) float32 786kB dask.array<chunksize=(4096,), meta=np.ndarray>
    lat       (cell) float32 786kB dask.array<chunksize=(4096,), meta=np.ndarray>
    lon       (cell) float32 786kB dask.array<chunksize=(4096,), meta=np.ndarray>
Attributes: (12/22)
    Conventions:                  CF-1.8
    atm_initial_conditions_file:  /lustre/orion/cli115/world-shared/e3sm/inpu...
    averaging_frequency:          1
    averaging_frequency_units:    nhours
    averaging_type:               INSTANT
    case:                         cess-control.ne1024pg2_ne1024pg2.F2010-SCRE...
    ...                           ...
    product:                      model-output
    realm:                        atmos
    run_t0:                       2019-08-01-00000
    source:                       E3SM Atmosphere Model (EAMxx)
    topography_file:              /lustre/orion/cli115/world-shared/e3sm/inpu...
    username:                     noel
scream_precip = scream.pr
scream_precip
scream_precip_mmhr = scream_precip * 36e5

plot 2D

import cftime 
start_time = cftime.DatetimeNoLeap(2020, 6, 1, 0, 0, 0, 0)
end_time = cftime.DatetimeNoLeap(2020, 6, 30, 0, 0, 0, 0)
scream_precip_june =  scream_precip_mmhr.sel(time=slice(start_time, end_time)).mean(dim="time").compute()
def get_NCL_colormap(cmap_name, extend='None'):
    '''Read an NCL RGB colormap file and convert it to a Matplotlib colormap
    object.

    Parameter:
        cmap_name     NCL RGB colormap name, e.g. 'ncl_default'
        extend        use NCL behavior of color handling for the colorbar 'under'
                      and 'over' colors. 'None' or 'ncl', default 'None'

    Description:
        Read the NCL colormap and convert it to a Matplotlib Colormap object.
        It checks if the colormap file is already available or use the
        appropriate URL.

        If NCL is installed the colormap will be searched in its colormaps
        folder $NCARG_ROOT/lib/ncarg/colormaps.

        Returns a Matplotlib Colormap object.
    '''
    from matplotlib.colors import ListedColormap
    import requests
    import errno
    import os, re
    import numpy as np

    #-- NCL colormaps URL
    NCARG_URL = 'https://www.ncl.ucar.edu/Document/Graphics/ColorTables/Files/'

    #-- read the NCL colormap RGB file
    colormap_file = cmap_name+'.rgb'
    cfile = os.path.split(colormap_file)[1]

    if os.path.isfile(colormap_file) == False:
        #-- if NCL is already installed
        if 'NCARG_ROOT' in os.environ:
            cpath = os.environ['NCARG_ROOT']+'/lib/ncarg/colormaps/'
            if os.path.isfile(cpath + cfile):
                colormap_file = cpath + cfile
                with open(colormap_file) as f:
                    lines = [re.sub('\s+',' ',l)  for l in f.read().splitlines() if not (l.startswith('#') or l.startswith('n'))]
        #-- use URL to read colormap
        elif not 'NCARG_ROOT' in os.environ:
            url_file = NCARG_URL+'/'+cmap_name+'.rgb'
            res = requests.head(url_file)
            if not res.status_code == 200:
                print(f'{cmap_name} does not exist!')
                raise FileNotFoundError(errno.ENOENT, os.strerror(errno.ENOENT), cmap_name)
            content = requests.get(url_file, stream=True).content
            lines = [re.sub('\s+',' ',l)  for l in content.decode('UTF-8').splitlines() if not (l.startswith('#') or l.startswith('n'))]
        else:
            raise FileNotFoundError(errno.ENOENT, os.strerror(errno.ENOENT), colormap_file)
    #-- use local colormap file
    else:
        with open(colormap_file) as f:
            lines = [re.sub('\s+',' ',l)  for l in f.read().splitlines() if not (l.startswith('#') or l.startswith('n'))]

    #-- skip all possible header lines
    tmp  = [l.split('#', 1)[0] for l in lines]

    tmp = [re.sub(r'\s+',' ', s) for s in tmp]
    tmp = [ x for x in tmp if ';' not in x ]
    tmp = [ x for x in tmp if x != '']

    #-- get the RGB values
    i = 0
    for l in tmp:
        new_array = np.array(l.split()).astype(float)
        if i == 0:
            color_list = new_array
        else:
            color_list = np.vstack((color_list, new_array))
        i += 1

    #-- make sure that the RGB values are within range 0 to 1
    if (color_list > 1.).any(): color_list = color_list / 255

    #-- add alpha-channel RGB -> RGBA
    alpha        = np.ones((color_list.shape[0],4))
    alpha[:,:-1] = color_list
    color_list   = alpha

    #-- convert to Colormap object
    if extend == 'ncl':
        cmap = ListedColormap(color_list[1:-1,:])
    else:
        cmap = ListedColormap(color_list)

    #-- define the under, over, and bad colors
    under = color_list[0,:]
    over  = color_list[-1,:]
    bad   = [0.5, 0.5, 0.5, 1.]
    cmap.set_extremes(under=color_list[0], bad=bad, over=color_list[-1])

    return cmap
<>:40: SyntaxWarning: invalid escape sequence '\s'
<>:49: SyntaxWarning: invalid escape sequence '\s'
<>:55: SyntaxWarning: invalid escape sequence '\s'
<>:40: SyntaxWarning: invalid escape sequence '\s'
<>:49: SyntaxWarning: invalid escape sequence '\s'
<>:55: SyntaxWarning: invalid escape sequence '\s'
/glade/derecho/scratch/weiyiwang/tmp/ipykernel_29936/248976764.py:40: SyntaxWarning: invalid escape sequence '\s'
  lines = [re.sub('\s+',' ',l)  for l in f.read().splitlines() if not (l.startswith('#') or l.startswith('n'))]
/glade/derecho/scratch/weiyiwang/tmp/ipykernel_29936/248976764.py:49: SyntaxWarning: invalid escape sequence '\s'
  lines = [re.sub('\s+',' ',l)  for l in content.decode('UTF-8').splitlines() if not (l.startswith('#') or l.startswith('n'))]
/glade/derecho/scratch/weiyiwang/tmp/ipykernel_29936/248976764.py:55: SyntaxWarning: invalid escape sequence '\s'
  lines = [re.sub('\s+',' ',l)  for l in f.read().splitlines() if not (l.startswith('#') or l.startswith('n'))]
cmp=get_NCL_colormap('WhiteBlueGreenYellowRed')

plot using imshow

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cf
import cmocean

# Example projection
projection = ccrs.Robinson(central_longitude=180)
fs = 24

# Create figure with 2 subplots
fig, axes = plt.subplots(
    ncols=1,
    figsize=(14, 5),
    subplot_kw={"projection": projection},
    constrained_layout=True
)

axes.set_global()
fg = egh.healpix_show(scream_precip_june, ax=axes, cmap=cmp, vmin=.1, vmax=.6)#cmocean.cm.thermal, vmax=0.5)
axes.add_feature(cf.COASTLINE, linewidth=0.8)
axes.add_feature(cf.BORDERS, linewidth=0.4)
axes.set_title("SCREAM Precip (zoom=7) June 2020", fontsize = fs )
fig.colorbar(fg)
plt.show()
../_images/37fd43fcf9743248a97f7acdb8aef04a6a533623da186eef91d7378e748e1fc1.png
# Create figure with 2 subplots
fig, axes = plt.subplots(
    ncols=2,
    figsize=(14, 5),
    subplot_kw={"projection": projection},
    constrained_layout=True
)

# Set same extent for both plots
extent = [90, 124, 20, 55]

# First subplot
axes[0].set_global()
axes[1].set_extent(extent, crs=ccrs.PlateCarree())
egh.healpix_show(scream_precip_june, ax=axes[0], cmap=cmp, vmin=.1, vmax=.6)
axes[0].add_feature(cf.COASTLINE, linewidth=0.8)
axes[0].add_feature(cf.BORDERS, linewidth=0.4)
#axes[0].set_title("SCREAM", fontsize = fs )
fg = egh.healpix_show(scream_precip_june, ax=axes[1], cmap=cmp, vmin=.1, vmax=.6)#cmocean.cm.thermal, vmax=0.5)
axes[1].add_feature(cf.COASTLINE, linewidth=0.8)
axes[1].add_feature(cf.BORDERS, linewidth=0.4)
fig.colorbar(fg)
plt.show()
../_images/712d2bc5c41f39145475503d3b118901d768fbe01e8621b3f31ee9d2fb2283bc.png

plot with contour

projection = ccrs.Robinson(central_longitude=180)
fs = 24

# Create figure with 2 subplots
fig, axes = plt.subplots(
    ncols=1,
    figsize=(14, 5),
    subplot_kw={"projection": projection},
    constrained_layout=True
)

extent = [90, 124, 20, 55]
axes.set_extent(extent, crs=ccrs.PlateCarree())
fg = egh.healpix_show(scream_precip_june, ax=axes, cmap=cmp, vmin=.1, vmax=.6, method="linear")
cn = egh.healpix_contour(scream_precip_june, ax=axes, levels=[0.3,0.6], colors="k", linewidths=1, alpha=1)
axes.add_feature(cf.BORDERS, linewidth=0.8)
axes.add_feature(cf.COASTLINE, linewidth=0.8)
axes.set_title("SCREAM Precip (zoom=7) June 2020", fontsize = fs )
fig.colorbar(fg)
plt.clabel(cn, inline=True, fontsize=8, colors="k", fmt="%0.1f")
plt.show()
../_images/c9466fb9da63453848c7c6192591335b69b92581c5a7fb6b08e849a9a8c5c8af.png

plot using UXarray

import uxarray as ux
uxds = ux.UxDataset.from_healpix(scream)
uxds
<xarray.UxDataset> Size: 15GB
Dimensions:   (crs: 0, time: 9457, n_face: 196608)
Coordinates:
  * crs       (crs) float64 0B 
  * time      (time) object 76kB 2019-08-01 00:00:00 ... 2020-09-01 00:00:00
Dimensions without coordinates: n_face
Data variables:
    pr        (time, n_face) float32 7GB dask.array<chunksize=(24, 196608), meta=np.ndarray>
    rlut      (time, n_face) float32 7GB dask.array<chunksize=(24, 196608), meta=np.ndarray>
    ELEV      (n_face) float32 786kB dask.array<chunksize=(4096,), meta=np.ndarray>
    LANDFRAC  (n_face) float32 786kB dask.array<chunksize=(4096,), meta=np.ndarray>
    lat       (n_face) float32 786kB dask.array<chunksize=(4096,), meta=np.ndarray>
    lon       (n_face) float32 786kB dask.array<chunksize=(4096,), meta=np.ndarray>
scream_precip = uxds.pr
scream_precip
scream_precip_mmhr = scream_precip * 36e5
start_time = cftime.DatetimeNoLeap(2020, 6, 1, 0, 0, 0, 0)
end_time = cftime.DatetimeNoLeap(2020, 6, 30, 0, 0, 0, 0)
scream_precip_june =  scream_precip_mmhr.sel(time=slice(start_time, end_time)).mean(dim="time").compute()
projection = ccrs.Robinson(central_longitude=180)

scream_precip_june.plot(
    projection=projection,
    cmap=cmp,
    features=["borders", "coastline"],
    title="SCREAM Precip (zoom=7) June 2020",
    width=700,
    clim=(0.1,0.6)
)
WARNING:param.GeoOverlayPlot01127: Due to internal constraints, when aspect and width/height is set, the bokeh backend uses those values as frame_width/frame_height instead. This ensures the aspect is respected, but means that the plot might be slightly larger than anticipated. Set the frame_width/frame_height explicitly to suppress this warning.
lon_bounds = (80, 124)
lat_bounds = (20, 55)

scream_precip_june_subset = scream_precip_june.subset.bounding_box(lon_bounds, lat_bounds)
#scream_precip_june_subset = scream_precip_mmhr.sel(time=slice(start_time, end_time)).subset.bounding_box(lon_bounds, lat_bounds).mean(dim="time").compute()
projection = ccrs.Robinson(central_longitude=180)

scream_precip_june_subset.plot(
    projection=projection,
    cmap=cmp,
    features=["borders", "coastline"],
    title="SCREAM Precip (zoom=7) June 2020",
    width=1100,
    clim=(0.1,0.6)
)
WARNING:param.GeoOverlayPlot03029: Due to internal constraints, when aspect and width/height is set, the bokeh backend uses those values as frame_width/frame_height instead. This ensures the aspect is respected, but means that the plot might be slightly larger than anticipated. Set the frame_width/frame_height explicitly to suppress this warning.